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We consider a quantum mechanical system represented in phase space (referred to 
hereafter as ’’Wigner space”), coupled to a harmonic oscillator bath. We derive quan¬ 
tum hierarchal Fokker-Planck (QHFP) equations not only in real time, but also in 
imaginary time, which represents an inverse temperature. This is an extension of a 
previous work, in which we studied a spin-boson system, to a Brownian system. It 
is shown that the QHFP in real time obtained from a correlated thermal equilibrium 
state of the total system possess the same form as those obtained from a factorized ini¬ 
tial state. A modified terminator for the hierarchal equations of motion is introduced 
to treat the non-Markovian case more efficiently. Using the imaginary-time QHFP, 
numerous thermodynamic quantities, including the free energy, entropy, internal en¬ 
ergy, heat capacity, and susceptibility can be evaluated for any potential. These 
equations allow us to treat non-Markovian, non-perturbative system-bath interac¬ 
tions at finite temperature. Through numerical integration of the real-time QHFP 
for a harmonic system, we obtain the equilibrium distributions, the auto-correlation 
function, and the first- and second-order response functions. These results are com¬ 
pared with analytically exact results for the same quantities. This provides a critical 
test of the formalism for a non-factorized thermal state, and elucidates the roles of 
fluctuation, dissipation, non-Markovian effects, and system-bath coherence. Employ¬ 
ing numerical solutions of the imaginary-time QHFP, we demonstrate the capability 
of this method to obtain thermodynamic quantities for any potential surface. It is 
shown that both types of QHFP equations can produce numerical results of any de¬ 
sired accuracy. The FORTRAN source codes that we developed, which allow for the 
treatment of Wigner space dynamics with any potential form, (TanimuranFP15 and 
ImTanimuranFP15) are provided as supplementary materials. 
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I. INTRODUCTION 


A Brownian oscillator (BO) model, which consists of a primary system coupled to a 
harmonic oscillator bath, is a versatile model that has been used to investigate fundamental 
problems in physics, chemistry and biology." The key feature of the Brownian model is 
that it describes irreversible dynamics through which the system evolves toward the thermal 
equilibrium state at finite temperature. This feature arises from interaction with the heat 
bath, which exhibits the canonical distribution at temperature T. To make the heat bath 
an unlimited heat source that possesses infinite heat capacity, the number of heat bath 
oscillators is effectively made infinitely large by replacing the spectral distribution of the 
system-oscillator coupling, J(cu), which was originally defined as the discretized distribution 
J{oj) = c^S(uj — ujj) (where Cj is the coupling strength between the system and the jth 
bath oscillator with frequency ujj), with a continuous distribution, for example, J[oj) oc c o. 
Because the time-evolution of the total system is described by the Schrodinger equation, 
the total energy is conserved and the dynamics are reversible. In the reduced description 
of the system obtained by tracing over the bath degrees of freedom using such methods as 
the path integral method^ or the projection operator method,- 4 '' however, the energy is no 
longer conserved, and its dynamics are irreversible, because the reduced system is merely a 
part of the total system. Heat bath effects arise in the reduced dynamics as fluctuation and 
dissipation in the reduced main system. These satisfy the classical or quantum version of 
the fluctuation-dissipation theorem. The reduced system evolves in an irreversible manner 
toward the thermal equilibrium state, in which the energy supplied by the fluctuations 
and the energy lost through dissipation are balanced, while the bath temperature does not 
change, because its heat capacity is infinite. 

With the above described features, the Brownian model exhibits wide applicability, de¬ 
spite its simplicity. This is because the influence of the environment can in many cases 
be approximated by a Gaussian process, due to the cumulative effect of the large num¬ 
ber of weak environmental interactions, in which case the ordinary central limit theorem 
is applicable,- while the distribution function of the harmonic oscillator bath itself also 
exhibits a Gaussian distribution. By adjusting the form of the spectral distribution, the 
properties of the bath can be adjusted to represent a variety of environments consisting 
of, for example, solid state materials, solvates, and protein molecules. This model has 
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been used to solve various problems of practical interest, in particular to investigate tun¬ 
neling pro cesses chemical reaction,—^ non-adiabatic transition, 13 ’ 14 quantum device 

systems,— ratchet rectification, 16,1 - to evaluate the efficiency of SQUID rings,— ^ and to 
analyze the line shapes in laser spectra.—^ 

While the Brownian model itself is fairly simple, it is somewhat difficult to apply in 
the quantum mechanical case not only analytically but also numerically, due to the infi¬ 
nite number of bath degrees of freedom. Analytically exact solutions of Green’s function 
for the BO Hamiltonian have been obtained only in the cases of a harmonic oscillator,a 
free particle,— and a free rotator,— using the path integral approach. Several approximate 
approaches have been developed to facilitate application of the Brownian model to more com¬ 
plicated systems. These approaches involve variational methods to study polarons^i— and 
the optical response of an anharmonic oscillator— using a damped oscillator as a trial func¬ 
tion, an instanton method for estimating the tunneling rate using instantaneously jumping 
paths between tunneling wells,a WKB method for evaluating the density matrix along a 
classical minimal action path,—— and diagrammatic expansion methods to study the anhar- 
monicity of potentials and the nonlinearity of the system-bath coupling.—The analytical 
expressions obtained in these studies are helpful to gain insight into the role of dissipative 
environments in the dynamics of systems, but they do not allow us to study situations in¬ 
vestigated in modern experiments that are usually described by complex potentials driven 
with time-dependent external forces. 

A great deal of effort has been dedicated to numerically calculating the time evolution 
of BO systems under external perturbations. Widely used approaches employ a reduced 
equation of motion that can be derived from the quantum Liouville equation with the full 
Hamiltonian by reducing the heat bath degrees of freedom. To obtain reduced equations 
of motion in a compact form, one usually employs the Markovian assumption, in which 
the correlation time is very short in comparison to the characteristic time of the system 
dynamics. In this case, the noise can be regarded as white. The quantum Langevin equation 
and the quantum Fokker-Planck equation have been derived with the projection operator 
method and the path integral method, for example.— 

In the classical case, the Langevin equation^ and the Fokker-Planck (or Kramers) 
equation 3 '” 18 have proved to be useful in the treatment of transport problems, and they have 
even been included in algorisms employed in molecular dynamics simulations. However, the 
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applicability of the quantum forms of these equations is very limited, because they can¬ 
not be derived in a quantum mechanical framework without severe approximations and/or 
assumptions. For example, in the treatment of the quantum Langevin equation expressed 
in operator form, it is generally assumed that the antisymmetric correlation function of 
the noise is very short (the Markovian assumption) and positive. A similar Markovian 
assumption has been used in the treatment of the quantum Fokker-Planck equation. But 
iu order for these assumptions to be valid, the heat bath must be at a sufficiently high 
temperature, in which case most of the important quantum dynamical effects play a minor 
role. This implys that the Markovian assumption is incompatible with obtaining a quantum 
mechanical description of dissipative dynamics at low temperature.— 

An Ohmic spectral distribution is generally assumed to realize Markovian noise. As we 
show in Appendix B, however, even if the dissipation process is Markovian, the fluctuation 
process may not be, because it must satisfy the fluctuation-dissipation theorem.- For this 
reason, if we apply the equation of motion under Markovian assumption to low temperature 
systems, then the positivity of the probability distributions of the reduced system cannot be 
maintained.— This is a fundamental limitation, known as the ’’positivity problem,” which 
is particularly significant for the quantum master equation. 41-18 If the system is not time 
dependent and if the system Hamiltonian and the system-bath interaction commute, the 
time-convolutionless (TCL) master equation becomes exact.— For the time-dependent 
case and/or non-commuting case, however, this master equation is valid only to second 
order with respect to the system-bath interaction, and the positivity condition is again 
broken. As a method to preserve positivity, the rotating wave approximation (RWA), which 
modifies the interaction between the system and the heat bath, has been applied in order 
to put the equation of motion in the Lindblad form. However, the RWA alters the thermal 
equilibrium state and the dynamics of the reduced system. These changes are particularly 
large in the case of a strong system-bath coupling and at low temperature. Moreover, in 
a typical quantum transport problem, the system is described by continuous energy states, 
and the energy levels of the heat bath and the system overlap. For this reason, the RWA 
cannot be used. Treatments of these kinds are therefore not sufficient to construct fully 
quantum mechanical descriptions of broad validity. 

Path integral Monte Carlo simulations do not have any of the limitations of the approaches 
discussed above, but this approach is computationally intensive, because the number of paths 
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to be evaluated grows rapidly with time, while sampling fails, due to the phase cancellation 
of wave functions.— - — Much effort has been made to overcome these problems and extend 
the applicability of this method.— - — Because this approach can easily incorporate the semi- 
classical approximation for the bath, it may be advantageous in the study of polyatomic 
systems treated in multi-dimensional coordinates, but applications to this point incorporat¬ 
ing full quantum mechanical dynamics have been limited to relatively small systems without 
time-dependent external force. 

Wave function based methodologies for the full Hamiltonian have been developed in order 
to avoid the reduced description of the system. The multi-configurational time-dependent 
Hartree (MCTDH) approach- 7 - employs time-dependent basis sets to represent the total 
wave function. Then, a variational principle is applied to derive the optimal equation of 
motion in order to reduce the bath degrees of freedom. This approach can be used to treat 
nonlinear system-bath coupling and anharmonic bath modes.— However, the number of bath 
modes must be increased until convergence is reached. This implies that the study of long 
time behavior requires more basis sets, which makes the calculation more difficult. In the 
effective-mode approach, the heat bath degrees of freedom are mapped to a linearly coupled 
harmonic oscillator chain. Then, the dynamics of the system are described by the wave func¬ 
tion of the system with a finite number of chained oscillators using a truncation scheme 13- ' 5 
or by utilizing the density matrix renormalization group method.— Strictly speaking, the 
time evolution obtained with the wave function based approach describes time-reversible 
processes and thus, within this approach, there exists no thermal equilibrium state. How¬ 
ever, in practice, this kind of approach has wider applicability than the reduced equation 
of motion. At this stage, the results obtained from these approaches have been limited to 
relatively simple systems. In particular, the inclusion of time dependent external forces is 
not as straightforward in these approaches as in the case of reduced equation of motion, 
because the energy of the total system changes due to the presence of an external force if 
the perturbation is strong, and hence the optimal basis set may also be changed. 

The reduced hierarchal equations of motion (HEOM), which are derived by differentiating 
the reduced density matrix elements defined by path integral, are reduced equations of mo¬ 
tion that can describe the dynamics of the system for non-perturbative and non-Markovian 
system-bath interactions with any desired accuracy under strong time-dependent perturba¬ 
tions at finite temperature.- In this formalism, the effects of higher-order non-Markovian 
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system-bath interactions are mapped into the hicrarchal elements of the reduced density 
matrix. In their original formulation, these equations of motion were limited to the case in 
which the spectral distribution function takes the Drude form (i.e., the Ohmic form with 
a Lorentzian cut-off) and the bath temperature is high." However, with the inclusion of 
low temperature corrections terms, this temperature limitation has been eliminated.— - — 
In addition, with the extension of the dimension of the hierarchy, in its present form, this 
approach is capable of treating a great variety of spectral distribution functions.— - — This 
formalism is valuable because it can be used to treat not only strong system-bath coupling 
but also quantum coherence between the system and bath, which is essential to study a sys¬ 
tem subject to a time-dependent external forced and nonlinear response functions.— - — The 
system-bath coherence becomes particularly important if the bath interaction is regarded 
as non-Markovian, as found from femtosecond nonlinear optical measurements, which are 
carried out on time scales that are much shorter than the noise correlation time of environ¬ 
mental molecules.— 

For a Brownian system, the reduced hierarchal equations of motion are expressed in the 
Wigner space representation.— - — In the Markovian limit, these equations of motion reduce 
to the Caldeira-Leggett quantum Fokker-Planck equation,—^ and in the classical limit, they 
reduce to the classical Fokker-Planck (Kramers) equation. 3 " 38 

Recently, the author derived the HEOM not only in real time, but also in imaginary 
time, which represents an inverse temperature, starting from correlated initial conditions for 
a system described by discretized energy states.— Reduction of these HEOM to a system 
represented in Wigner space is not straightforward, because they involve derivatives with 
respect to the position and momentum that require a careful treatment with regards to the 
order of time slices in the path integral formalism. In this paper, we present the derivation 
of real- and imaginary-time HEOM in Wigner space and demonstrate the validity of these 
equations. 

The organization of the paper is as follows. In Sec. II, we present a model Hamiltonian 
and its influence functional with correlated initial conditions. In Sec. Ill, we derive the real¬ 
time quantum hierarchal Fokker-Planck (the real-time QHFP) equations using the influence 
functional given in Sec. II. In Sec. IV, we derive the imaginary-time quantum hierarchal 
Fokker-Planck (the imaginary-time QHFP) equations, which are convenient for evaluating 
thermodynamic quantities of the system. In Sec. V, the validity of our approach is demon- 
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strated through numerical integration of the real- and imaginary-time QHFP equations for a 
harmonic system and comparing the calculated results with the exact results obtained from 
analytical calculations. Section VI is devoted to concluding remarks. 

II. REDUCED HIERACHAL EQUATIONS OF MOTION FROM 
CORRELATED INITIAL CONDITIONS 

We consider the situation in which the system interacts with a heat bath that gives rise 
to dissipation and fluctuation in the system. To illustrate this, let us consider a Brownian 
Hamiltonian expressed asA- 

Htot = H A (p, q) + 

j 

where 

^2 

HA(j>,q) = ^ + m ( 2 ) 

is the Hamiltonian for the system with mass m and potential U ( q ) described by the momen¬ 
tum p and position q. The bath degrees of freedom are treated as an ensemble of harmonic 
oscillators, and the momentum, position, mass, and frequency of the jth bath oscillator are 
given by pj, Xj, rrij and Uj, respectively. In the conventional Brownian model, the system- 
bath interaction is represented by a bilinear function of the system and bath coordinates as 
Hi = -«E, oijXj. Brownian models employing this bilinear interaction have been studied 
with various approaches.— In this paper also we restrict our investigation to this bilin¬ 
ear form to simplify the derivation, but we note that extension to the non-bilinear case is 
possible.- 1 ^— To maintain translational symmetry in the case of U(q) = 0, required to de¬ 
scribe the motion of a free Brownian particle, we include the counter-term JU a^q 2 /2m,ji0j 
in Eq.(H]). 

The heat bath we consider is characterized by the spectral distribution function defined by 
= J2j(ha‘j/2mjUjj)5(uj — coj) and the inverse temperature, /3 = l/ksT, where ks is the 
Boltzmann constant. The path integral used here to derive the HEOM is expressed in terms 
of an influence functional with correlated initial conditions. The influence functional that we 
employ, Fci[t,f3h], is calculated by taking the trace over the heat bath degrees of freedom, 
starting from the thermal equilibrium state of the total Hamiltonian. The calculation of the 


Pj , 

2 rrij 2 


Xj 


ajq 

mjuj] 


( 1 ) 
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influence functional for a heat bath consisting of harmonic oscillators is analogous to that of 
the generating functional for a Brownian oscillator system if we regard the system operator 
in the system-bath interaction q as an external force acting on the bath.—As shown 
in Appendix A, the reduced density matrix elements of the system with correlated initial 
conditions can be expressed as 

i rq=q{t) rq'=q'(t ) rq' 0 =q(P h ) 

P{q,q']t)=— D[q(t)\ D[q'(t)\ D[q(r)} 

^tot jq 0 =q( 0) Jq' Q =q'(0) Jq 0 =q(0) 

x e i SA ^F CI [q, q', q; t, Ph]p eq [q ; ^h]e~^ SA (3) 


where Sa[ q] £] is the action for the Hamiltonian of the system, Eq. (l2]l . given by 


S A [q-,t]= dr < —mq 2 {r) — U{q{r)) > , 


(4) 


and p eq [q ; /3h] is the initial thermal equilibrium state, with the heat bath defined by Eq. (lA12lh 
We assume that the spectral density J(oj) has an Ohmic form with a Lorentzian cut-off 
and write^ 

hm( y 2 c o 


J M = 


(5) 


7T 7 2 + LU 2 ’ 

where the constant 7 represents the width of the spectral distribution of the collective 
bath modes and is the reciprocal of the correlation time of the noise induced by the bath. 
The parameter ( is the system-bath coupling strength, which represents the magnitude of 
damping. This spectral distribution approaches the Ohmic distribution, J(oj) ~ hm^uj/n, 
for large 7 . In Appendix B, we present several profiles of fluctuation and dissipation terms 
for the Drude distribution to illustrate the origin of the positivity problem in the Markovian 
master and Redfield equations. 

With J(oj) given by Eq. (]5]1 . the influence functional with correlated initial conditions is 
expressed as^ 

Fci[q, q', q; t, f)h}= <T Sl ’“"W *** ^(O+CoW-tS^)} 

-fqdt" E 

xe k=1 , ( 6 ) 

where, for the Matsubara frequency Uk = 2 tt k//3h, we have defined 


®(t) = ft lv(t) - q'(t)}, 


( 7 ) 






( 8 ) 


©o(*) = | [g(t) + g'(*)] - ^7 cot 0-^- ) [q(t) - q\t )] |> , 


Go( 0 ) = ^ [g( 0 ) + g'( 0 )] , 


(9) 


P 


r ph f i oo 

l iTW) {^% 


e (W . 2 ^ 2 7"—/ 1 


[7 cos(z/ fc r / ) - ii/ fc sin(i/ fc r')] [ 
7 2 - *2 J ’ 


and for k > 1 , 


e t (0 ^ (,(() - g'W!, 


( 10 ) 


(ii) 


r<jr' g V) ^ |c0s( ^ ) - i 2 sill(l,tT,)1 . 
P Jo t - K 


and 




E 


2 7 2 


T-*'* 


2^k 


W) - q'{t )] 2 , 


( 12 ) 


(13) 


L k=K+l 

where + cu 2 ) is the correction factor that counteracts the overestimation of the 

contribution of higher-order Matsubara frequencies approximated by the delta function with 
cut-off number, K , introduced in Appendix B for the characteristic frequency of the system, 
oj c . This modification improves the convergence of hierarchies at lower temperature. We 
now introduce the hierarchal elements that play an essential role in our formalism: 

1 


A 


q=q(t) rq'=q'(t) rq'o=q(Ph) 

Dm ] / D[q'(t)} / D\q{r)\ 

tot Jq 0 =q( 0 ) Jq'o=q'{ 0 ) ^<? 0 =q( 0 ) 


x^ SA[q ’ t] F^..., jK [q, q', g; t, mp eq W, m*-* SAWA , 


(14) 


where 


A; t,^]= 


- J e -7 * 


[yo 


dt'e 7 *' 760(0 + G o( 0 ) - -0(/3/i) 


x 


n 

fc=i 




bo 


dfVA fc 0 fc (f') - yT fc (/3£) 
n 


Jk 


*-F c i[q , 5 ',?; 


(15) 


for nonnegative integers n,j 1 ,..., Jjk- From the above definition, the first hierarchal element 
and the reduced density matrix given by Eq.(j3]) are identical: p(q,q']t) = Po^..o(q,q l ’,t). As 
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shown in Appendix C, we then have the following equations of motion: 


d Phi,3K M’t) 


dt 


K 


-L (q, q') + ny + ^j k Vk + E'(g, q') 


k =1 


PnljK b’M 


w70o(g,? , )Pi"... 1 ] ir fog';*)) 

K 


J^J^kQkiq, q )p-i-h i. ,/k («> e'; *) 

fc=i 

*(«> «') (s> <?'; *) + S (s> 9'; *)), 


k =1 


where 


£ (9, 9') = 


©o(g,9) = — 


h 2 d 2 h 2 <9 2 


2m dq 2 

2m dq' 2 

( 9 - 

d 

\ | my 

\dq 

dq' 

) h 


+ C/(9)-C/(9'), 
/^7 


cot 


(9 -«') 


(16) 


(17) 


(18) 


and <f>(g, </), @fc(g, g 7 ), and S'(g,g') are defined by Eqs. (171). (fill) , and (fT3j) by making the 
replacements q(t) —y q and q'{t) —y q'. In the HEOM formalism, only the first element 
p(q,q']t) = p { ol,o(q,q']t) h as a physical meaning and the other elements ( 9 , q'\ t) 

are introduced in numerical calculations in order to treat the non-perturbative and 11011 - 
Markovian system-bath interaction. We can evaluate p^ 0 (q, q';t) through numerical inte¬ 
gration of the above equations. 

We next explain the truncation scheme that we use for the hierarchical equations, which 
is different from the scheme used in previous studies.— - — First, we choose the number of 
Matsubara frequencies to be included in the HEOM, A", such that it satisfies K u)q/v\. 
Then, we introduce the scaled integer A " 7 as K 7 = mt(Kv\/'y) for zy > 7 and K 1 = K 
for U\ < 7 , which allows us to make calculations in the highly non-Markovian case more 
efficiently. The index for the hierarchy, denoted by n, for a given value of 7 , then runs from 
0 to K 7 . The total number of hierarchy members to be included in the calculations is then 
given by N = (AT, + K + 1)!/(AT + 1)!/A' 7 !. For the case ^2 k=1 j k > AT, we truncate the 
hierarchal equations by replacing Eq. ffTBl) with 


S n ) 


L + -) ( 19 ) 


In practice, we can simply set pi") Jk (q, q’\ t) — 0 instead of employing the above equation, 
because p(© JK (q, q'\ t) decays to zero as t becomes large.— For the K 1 and K 7 +1 members 
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of the hierarchy, we have the following relation, valid to order 5t: 


o 1} (9, 1 { -7©o(?, 9 ')pC-,o(9, 9'; *) 


/\ JW) 


iV+ 1 


$( 9 , 9 ') 


- -©o(9,9>S(9,9';i) 


Po A .T,o 2) (9, 9'; *) - 6!..(9, 9'; *) 


ic 


fe=i 


( 20 ) 


This asymptotic relation allows us to obtain the terminator for given 7 in the form 


dpo K ?l ( 9 , 9 '; t) 


dt 


h 


L (q, q') + Ayr - $(q, q')®o{q, 9 ') + -( 9 , 9 ) 


p£o (9, 9'; *) 


- A' 7 7 © 0 (g, 9 ')Po A T,o 1} (9, 9'; *) 


( 21 ) 


This equation reduces to the quantum Fokker-Planck equation in the Markovian limit, i.e., 
the Ohmic distribution (7 —> 00 ) with the high-temperature limit.—— 

While the terms 0 and from the correlated initial state do not appear in Eqs. ffl 6 )l 
and CD, they define the hierarchal elements for the correlated initial equilibrium state.— 
To demonstrate this point, we consider the initial states of the density operators, obtained 
by setting t — 0 in Eq. ffITj) : 


filjiM-'O) = ± (") (G,«r^(^;0), 

171=0 ' 


( 22 ) 


where 




l ri'o=<i(PK) 


Qo = q(o) 


/I- \ Jk 

mr)} ( - S ©(/3R)) n(- S *^»)) p[9l/?fi] 


fc=l 


h 


(23) 


are the equilibrium hierarchal elements. Here, Za, Z tot , and Zb are the partition functions 
of the system, total system, and bath, respectively, related as Za = Z tot /ZB■ We then have 
p eq [q-,fih} = Z B p[q ; /3h]. It is important to note that the steady state of j K (t) for n > 0 
in Eq. (ll 6 |) is slightly shifted from the initial thermal equilibrium state as a result of the 
influence of the sum in Eq. (1221) . However, because Po°?..o(0 is n °t influenced by this effect, 
expectation values calculated using p^ 0 (t) does not change. 

From the above definition, it is clear that the HEOM members at time t — 0 represent 
a correlated initial state, while the zeroth member, pfp 0 ( 0 ) = p[q;/3H], involves the static 
correlations. In Fig. [9l the correlations responsible for the correlated initial state are 
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represented by green arcs, and the static correlations are represented by red arcs. After the 
time evolution, the elements p^' j K (q,q'',t) describe the dynamical correlation, represented 
by the blue arcs and lines in Fig. [HI 


III. REAL-TIME QUANTUM HIERARCHAL FOKKER-PLANCK 
EQUATIONS 


We now introduce the Wigner distribution function, which is the quantum analog of the 
classical distribution function in phase space. For the density matrix element p'^[ j K (q, q'] t), 
this is defined as^i^r— 


W, 


(n) 


Jl.-JK 


(P, 91 t) 



X X 

2'“ + 2 ]t 


(24) 


The Wigner representation of the reduced density matrix defined in Eq. (l3]l . W(p,q;t), and 
the hrst member of the hierarchal elements are then identical: W(p,q\t) = Wq°_ 0 (p, q:t). 
The Wigner distribution function is a real function, in contrast to the complex density 
matrix. In terms of the Wigner distribution, the quantum Liouvillian takes the form—hi 


^QmKLjAP’ 9 ) = 


.L^w {n) ■ (p q)-- 
mdq 31 ....jkWV h 


l^ u w{p-p\ qWj.L.jJr- 9), 

(25) 


where Uw(p, q) is given by 

Uw (p, 1) = 2 f dx An (i [u (q + §) - U (9 - f ) 


The quantum Liouvillian can also be expressed a& 


113.114 


- q) = 


p d 1 f / h d \ / h d 

mdq^ ih\ \ 2 i dp) \ + 2 i dp 


(26) 


W, 


(n) 


n,—,jK 


ip, q)- 

(27) 


While the above expression is easier to integrate in the case that the potential is nearly 
harmonic, the expression in Eq. (J25]) is numerically stable, and it can be applied with any 
form of potential, including an unbounded potential. 

Using the Wigner distribution and quantum Liouvillian, the equations of motion appear¬ 
ing in Eq. (TT6|) can be expressed in the form of quantum hierarchal Fokker-Planck (QHFP) 
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equations in real time as 


d_ 

dt 


W. 


(n) 


3 1.— ,3K 


C V , q]t) = - 


+ $ 


A' 


£qm + ™7 + E JfcZ/fc + 


fc=l 


W. 


(n) 


31,— ,3k 


(p, A 




A 


(n) 


fc=l 


n'yO 0 W^ jK (p,q- 1 t) 


K 




jk— l,-” Pif 




(28) 


fc=l 


where <h = d/dp, 


@o = C 


mfoy 

p H-cot 


d ' 
dp 


and 


= 

rnQ 

~T 


2m / y 2 ( d 

2 7 2 


E 

\_k=K+l 


CV 


7 


z'n 


9p 2 


(29) 

(30) 

(31) 


As in the case of the energy eigenstate representation,— the above equations are identical to 
the equations derived from factorized initial conditions.— - — The above equations are then 
truncated by using the modified ’’terminators” expressed in the Wigner representation. As 
explained in Sec. II, the number of Matsubara frequencies to be included in the calculation, 
K, is chosen to satisfy K uj c / u\. The upper limit for the number of hierarchy members 
for given 7 is then chosen to be A,, = /y) for 17 > 7 and K 1 = K for 17 < 7 . Then, 

for the case Yjk=\ jk > A, we truncate the hierarchal equations by replacing Eq. (1281) with 

d rrr (n) „.+\ _ ( r 1 0A ta/ ( n ) -A ^32) 


foWju-jiciP’ 9; A = - [£qm + S'J W n ..... jR (p. q; t ), 
while, for the case n = A 7 we employ 


d 


w-<E(p, 9 ; A = - t-’QM + a 7 7 - $00 + s 7 MAA'oAh 9 ; A - a 7 7(-) 0 it 0 A;: 0 ,j (a 9 ; A 


(A)/ 


7a 7 -i) ( 


dt 


(33) 


We can evaluate ?/f (p, g; i) through numerical integration of the above equations. While 
only the first element W(p, g; t) = Wo!o,- o(P, 9 ; t) has a physical meaning and the other ele¬ 
ments W-/ } j k (p, g; t ) are initially introduced to avoid the explicit treatment of the inherent 


13 
















memory effects, it turns out, however, that these elements allow us to take into account the 
system-bath coherence,- ent anglement ffi 116 1 117 an d expectation values that include the bath 
operators as (Hi) = — (< 7 ^ 07 %)-“ The HEOM consist of an infinite number of equations, 
but they can be evaluated with the desired accuracy by depicting the asymptotic behavior 
of the hierarchal elements for different K and using this to determine whether or not there 
are sufficiently many members in the hierarchy. Essentially, the error introduced by the 
truncation to be negligibly small when K is sufficiently large. 

The correlated initial equilibrium state defined by Eq. (j22il is expressed in the Wigner 
representation accordingly. The correlated initial equilibrium state can be set in the HEOM 
formalism by running the HEOM program until all of the hierarchy elements reach the steady 
state and then use these elements as the initial state,- or by integrating the imaginary-time 
HEOM that we discuss in the next section.— In practice, the former approach is simpler, 
because it requires the real-time HEOM only. This approach has been used to set the 
correlated initial conditions of the HEOM derived from factorized initial conditions that are 
identical to those used with the present HEOM. 

The HEOM in Wigner space is ideal for studying quantum transport systems, because it 
allows the treatment of continuous systems, utilizing open boundary conditions and periodic 
boundary conditions .—' 101 In addition, the formalism can accommodate the inclusion of an 
arbitrary time-dependent external field . 95-97,105 

In the Markovian limit, 7 —>■ 00 , which is taken after the high temperature limit, yielding 
the condition /3h r 7 <C 1, we have the quantum Fokker-Planck equation^— 

9 ; t) = -t QM W^(p, q; t) +(^(p + ^ ( 0 | (p, 91 t), (34) 

which is identical to the quantum master equation without the RWA.- Because we assume 
that the relation <C 1 is maintained while taking the limit 7 —> 00 , this equation 
cannot be applied to low-temperature systems, in which quantum effects play a major role. 
As in the case of the master equation without the RWA, the positivity of the population 
distribution, P(q) = f dpW(p,q;t), cannot be maintained if we apply this equation in the 
low temperature case. 

The classical HEOM can be derived by taking h — y 0 .— The Wigner distribution 
function reduces to the classical one in this limit. The classical equation of motion is helpful, 
because knowing the classical limit allows us to identify the purely quantum mechanical 
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effects.— ' i 102 i 105 


IV. IMAGINARY-TIME QUANTUM HIERARCHAL FOKKER-PLANCK 
EQUATIONS 


The equilibrium reduced density matrix has bee n eva luated with several approaches/ 


118019 


By applying the methodology developed in Ref. 


108 


we can derive the quantum liierar- 


chal Fokker-Planck (QHFP) equations in imaginary time. This allows us to calculate the 
thermal equilibrium distribution W eq (p, q) at inverse temperature (3h. Instead of the quan¬ 
tum Liouvillian, this equation involves the left-sided operators Ha(P, q) and q. While the 
Wigner transformations of these operators become complex operators, the Wigner distri¬ 
bution W eq (p, q) is a real function. In order to make the numerical calculations easier to 
carry out, we rewrite Ap as {Ap + pA)/2 to perform the Wigner transformation, where A 
is an arbitrary operator. Other than this, the derivation of the imaginary-time QHFP is 
parallel to that of the imaginary-time HEOM in the energy eigenstate representation.— By 
introducing the Winger distribution for imaginary time, IF^r \ m (p, g; r), which is defined 
by the Wigner transformation of the density operator that in path integral form is given by 






m—l 

t )\ n 

J 9(0) = 90 g= 1 

2m 

* n 

g'=m—l +1 


dTgCOS(u k9 Tg)q(Tg 




'0 


Sin(i vvMv) I p[«. o'; t], 


(35) 


we obtain the imaginary-time QHFP equations as 

i K 

[m:l] 


dr 1 


( t )= -HaW^^t) + t Cfcm+1 ^s(u km+1 r)qW^i +1 (r) 


fcl fcm \ 


h 


fc m+l =0 


K 


c k m +1 sm(u k m.+ir)qW^ n ^ k 

km+l= 0 
m—l 


[m+l:Z+l] / \ 
m +1 V 1 ) 


^ cos(u k hr)qW^J k Li 


k h + 1 ,...,k rn 




h=1 


- yy sin {v k hT)qW } 


[m—l:l—l] 

fcl fch—l fch-\-l fcm 




(36) 


h=m-l -\-1 


where the factors c k are expressed as c 0 = mCj//3 and c k = 2m('y 2 //3(y + v k ), for k > 1. We 


set W, 


[mil] 


k 1 ,...,k rn 


Y) = 0 for higher-order elements in hierarchy denoted by m to truncate. The 
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Euclidean Liouvillian is expressed as 


R * W = i £ ~ j£) + 5 £ ^ »>• (37) 


with 


U'(p,q)=f dxsin(jp) |t/' (g + |) + U' (? ~ £) } > ( 38 ) 

for the potential, I/'(g) = I/(g) + m(^ n /q 2 /2, including the counter-term. This can also be 
expressed in differential form as 


H a 


A ( 2 _ 

2 m 4 d 2 q) 


1 

2 



/ h d 
\ q ~ Jidp 


+ u' 



(39) 


If the anharmonicity of the potential is small, the above expression is useful. The initial 
conditions p[° : °] (g, q) — 1 and pl 0; °] (g, g') = 0 for q ^ q' are expressed as W^(p : g; 0) = 1/27r. 
By integrating Eq. (136|) from r = 0 to r = f3h , we can evaluate the equilibrium distribution 
function W ( p , g; /3h). 

Once we obtain the equilibrium distribution, we can calculate the partition function 
employing the relation 


Z A {Ph) = I dp j dqW^°\p, q] ph). (40) 

This allows us to calculate the Helmholtz free energy, F A = — \ii(Za)/P, the entropy, 
S'a = ksfPdF a! dfi , the internal energy, U A = —9In(Z^)/9/3, and the heat capacity, 
Ca = —kB(3 2 dUA/d(3 for any potential. If the system is subject to an external force A f(p, g), 
where f(p, q) is any function of the momentum and position, p and g, we can also calculate 
the susceptibility, xa = —(dF/d A), from Z A . 

It should be noted that even if the potential is a function of time, we can calculate 
thermodynamic quantities as functions of time through Za(/3H ; t), assuming that the system 
reaches the thermal equilibrium state faster than the change of the potential. 


V. NUMERICAL RESULTS 

In principle, the HEOM provide an asymptotic approach that allows us to calculate 
various physical quantities with any desired accuracy by adjusting the number of hierarchal 
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elements. Here, we demonstrate the applicability and validity of the real-time and imaginary- 
time QHFP equations, by presenting the results obtained from numerical integrations of 
Eqs.fl2ED-dMD and Eqs. (13(j]) - (138|) . For this purpose, we consider the harmonic potential 



(41) 


From our numerical solutions of Eqs. (128|) - (j33]L we have computed the equilibrium distri¬ 
butions, the auto-correlation functions, the first- and second-order response functions and 
examined the roles of a non-factorized thermal state, and the roles of fluctuation, dissipation, 
and system-bath coherence. From those of Eqs. (I36l) - (l38l) . we have computed the equilib¬ 
rium distributions and thermodynamic quantities. Below, we compare these results with the 
same quantities calculated from analytically exact expressions for the Brownian oscillator 
system^ - and from the time-convolutionless (TCL) Redheld equation both with and without 
the rotating wave approximation (RWA)^r— (see Appendix D) as critical non-perturbative 
and non-Markovian tests. Note that the TCL equation is exact if the system Hamiltonian is 
time independent and if the system Hamiltonian and the system-bath interaction commute. 
However, here we consider the non-commuting case. 

Below we also present our results for calculations of thermodynamic quantities obtained 
from the imaginary-time QHFP and compared them with analytical results. 

A. Steady state distribution: Static system-bath coherence and mixed state 

For a harmonic system, the equilibrium distribution in the Wigner representation is 
analytically expressed as 



where N = 27Tis the normalization factor and ( q 2 ) and ( p 2 ) are the mean squares 


of the position and momentum, respectively. 

The Wigner distribution for an isolated oscillator is written W^(p,q). For the Hamilto¬ 
nian Eq. flUJ) , we have^ 



(43) 


and 



( 44 ) 
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(a) Quantum Hierarchal FP (b) TCL Redfield 



FIG. 1. (a) The initial conditions (blue curves) and steady state solutions (red curves) for the 
low temperature case j3h = 3.0, calculated from (a) the real-time QHFP, (b) the TCL Redfield 
equation, and (c) the TCL Redfield equation with the RWA. The other parameter values are 
wo = 1.0, 7 = 1.0, and ( = 1.0. The factorized initial state given by W^(p, q) with Eqs. (l43l) and 
(1441) is set as the temporally initial state at time t = 0. After integrating the real-time QHFP 
and the TCL Redfield equations for a sufficiently long time (t = 100), the distribution reaches the 
steady state. In the real-time QHFP case, the obtained steady state is identical within numerical 
error to the thermal equilibrium state W^ 0 {p,q) with ( q 2 ) and (p 2 ) given by Eqs. (H5jl and (l46j) . 
while those from the TCL Redfield equations are similar to the original factorized initial state. 
This implies that the TCL Redfield equation cannot take into account the system-bath correlation 
properly. 


The Wigner distribution for a harmonic Brownian system is denoted by W^ g 0 (p,q). In this 
case, we have^- 

1 A 1 


(q )bo — 


m P 77 7> + "I + l‘5r 2 K)| 


(45) 


and 


<P 2 >bo = £ £ 


u,i + \ir 2 ( Vk )\ 


P + \ sr2 M\’ 


(46) 


with AT 2 (u;) = ( n f 2 uj/('y 2 + ui 2 ). In the Wigner representation, the thermal equilibrium state 
under the factorized assumption, exp(— /3Ha) exp(— (3Hb), is denoted by W^(p, q ), while the 
true thermal equilibrium state of the reduced density operator, tr B {exp[— /3(Ha+Hj+H b )}}, 
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is denoted by W^ q 0 (p,q); the difference between the two distributions arises from the static 
system-bath coherence and represents the non-factorized effect of the thermal equilibrium 
state. 

To obtain the thermal equilibrium state from the real-time QHFP, we integrated Eqs. 
m - m f rom a temporal initial state until all of the hierarchy elements reach the steady 
state. In principle, the initial state can have any form, but to elucidate the difference 
between the factorized (pure) equilibrium state and the true correlated (mixed) equilibrium 
state, we chose 0 (p, g; 0) = W A (p,q) and ]K (p, q\ 0) = 0 for other elements in 

the QHFP case. For the TCL Redfield case, we chose pjj( 0) = exp (—f3E'-)/Z' A , where 
Ej is the jth eigenenergy of Eq. (l4TT) with the counter-term H' A = ( Ha + m(^q 2 /2), and 
Z' A = exp(— (3Ej), as explained in Appendix D. 

For all of our computations, we fixed the oscillator frequency as ujq = 1.0. Then, we chose 
the coupling strength, inverse correlation time, and inverse temperature as ( = 1 , 7 = 1 , and 
(3 = 3. We thus consider the case of intermediate coupling strength and low temperature. 
For the QHFP, we set K = 7, which leads to the depth in terms of 7 as AQ = 2 and the 
total number of hierarchy elements N = 4268. The mesh size of the Winger function was 
optimized for the Liouvillian given in Eq. (l25lh — and we used n q = 80 and n p = 30 for 
the region \q\ < 2.8 . For the TCL Redfield equation, we employed six eigenstates. The 
calculated results and factorized initial state were translated into the Wigner distribution 
through Eqs. (1D8|) and (ID9D . respectively. 

In Fig{l](a), we display W { f y [ 0 (p, q\ t) for the factorized initial state at t — 0 given by 
W A (p,q) (blue curves) and the steady state distribution at t = 100 (red curves) obtained 
from the real-time QHFP calculation. We found that even if we start from the factor¬ 
ized initial state, the steady state solution is the true thermal equilibrium state, denoted 
by W^ q 0 (p, q). This indicates that the real-time QHFP has the capability to produce the 
thermal equilibrium state with a static system-bath correlation through the fluctuation and 
dissipation terms. In the TCL Redfield equation cases, Figs. mb) and (c), the calculated 
steady states (red curves) are similar to the factorized initial states (blue curves). The peak 
intensity of the TCL result in the case without the RWA is slightly higher than that in the 
case with the RWA, because the ground and first excited populations in the former case 
are poo(t ) = 1.083 and puit) = —0.090, due to the breakdown of the positivity condition, 
which is a physical requirement for the reduced equations of motion necessary for the pop- 
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illation state of the density matrix to be positive. 41-48 Other than this difference, the TCL 
distributions are similar to the factorized distribution. This indicates that the TCL Red- 
field equation cannot take into account the static system-bath coherence, because the TCL 
theory in the present case is valid only to second order in the system-bath coupling. 

As explained in Appendix A, the effects of the system-bath coherence consist of the 
imaginary-time (static) part and complex-time (correlated) part represented by the red and 
green arcs in Fig. [9l respectively. From the equilibrium distribution, we can only observe 
the effects of the static part. To elucidate the correlated part, we need to calculate the 
nonlinear response function, as will be discussed in Sec. V-C. 


B. Two-body correlation functions: The roles of fluctuation, dissipation, and 
non-Markovian effects 


We next calculate the two-body correlation functions to investigate the roles of dissi¬ 
pation, fluctuation and non-Markovian dynamics. The symmetric correlation and linear 
(first-order) response functions of the position are defined by C{t) = ( q{t)q + qq(t))/2 and 
i?,d)(t) = ([q(t), q\)/h, respectively. While the auto-correlation function of the position is 
given by C(t), the observable of a linear measurement involving infrared and THz spectra, 
which are expressed in terms of a dipole proportional to q corresponds to R ^ ( t ). The Fourier 
transformation of these functions are denoted by C[u\ and They are expressed as the 

real and imaginary parts of the normalized spectral distribution for the Brownian oscillator 


J'(co) = 


h 


m {oJq — u> 2 ) + i(jl[i(jj] ’ 
where I[s] is the Laplace transformation of B(t) defined by Eq. (1A8D as 


r 00 1 

I[s] — / dt—B(t)exjp(—st). 

Jo 


m 


For the Drude distribution, Eq.Q, we have I[s] = C7/( s + 7 ) an d 

h 5r 2 (o;)coth(^) 


CM = 


m (a ,2 _ a;2 _ + (5r 2 (a;))^ 


and 


= 


h 


5r 2 (o;) 


m (cu 2 - a; 2 - 5fi 2 (a;))' + (<W 2 (a;)) 


2 > 


(47) 


(48) 


(49) 


(50) 
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where 5n 2 (u) = Cl^Kl 2 + w 2 ). 

In order to calculate the above functions using an equation of motion approach, we employ 
the following forms:— 

c (t) = \ tr {qG(t)q°pz} ( 51 ) 

and 

R { 1 ) (t) = j^tr | qG{t)q x pZ } , (52) 

where q x A = qA — Aq,q°A = qA + Aq, G(t)A = e - lH t°t t / h A e lHtott / n for any operator A, and 
Pm = e- pHtot /Z to t with Z tot = tr{pZ). 

In the reduced equation of motion approach, the density matrix is replaced by a reduced 
one. In the QHFP case, pZ is replaced by the hierarchy member (p,q;t), whereas 

in the TCL Redheld case, it is replaced by pjk(t). The Liouvillian in G{t) is replaced using 
Eqs. Q2BD -QMD and Eqs. (1D2D - (ID4I) . respectively. 

We evaluate Eqs. (1521) and (15T]) in the following Eve steps.- 1 ^ (i) We hrst run the 
computational program to evaluate Eqs. (I28l) - (l33l) in the QHFP case and Eqs. (1D2|) - (1D4I) 
in the TCL Redheld case for sufficiently long times from the temporal initial conditions to 
obtain a true thermal equilibrium state, as illustrated in Sec. V-A. In the QHFP case, the 
full hierarchy members j K (p,q', 0 ) are then used to set the correlated initial thermal 

equilibrium state, (ii) The system is excited by the hrst interaction q x or q° at t = 0. In the 
Wigner representation, they are expressed as d/dp and 2 q, respectively, (iii) The evolution 
of the perturbed elements is then computed by running the program for the QHFP or TCL 
up to time t. (iv) Finally, the functions defined in Eqs. (l52|) and (I5T)) are calculated from the 
expectation value of q. By performing a fast Fourier transform, we obtain their spectra. 

In computing the results reported below, we chose the number of Matsubara frequencies 
for the QHFP equation as K = 5-8, which leads to the depth in terms of 7 as /C, = 3-6 
and the total number of hierarchy member N = 601-16093. The mesh size of the Winger 
function was optimized for the Liouvillian given by Eq. (125|) . n - and we used n q = 80-120 
for the region |g| < 4-6 and n p = 30-120 for the region |p| < 2.8-11.2. In the TCL Redheld 
cases with and without the RWA, we varied the number of energy levels between 6 and 16 
depending on the temperature. 
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FIG. 2. The auto-correlation (symmetric correlation) function of the Brownian oscillator system 
for several inverse temperatures: (a) j3h = 3.0, (b) /3h = 1.0, (c) f3h = 0.5. The dotted, red, blue, 
and blue-dash curves represent the results obtained from the analytic expression Eq. (|49l) . the 
QHFP, the TCL-Redfield, and TCL-Redfield with the RWA, respectively. The intensity of each 
line is normalized with respect to its maximum peak intensity. The other parameters values are 
fixed as ojq = 1 . 0 , 7 = 1.0 and £ = 1 . 


1. Auto-correlation function: Fluctuation and temperature effects 


First we study the temperature dependence of the auto-correlation function for the fixed 
coupling strength £ = 1 and the inverse noise correlation time 7 = 1. In Fig. [2j we 
compare the calculated real-time QHFP results obtained from Eqs. (j28]) - (l33|) with analytical 
results obtained from Eq. (l49|) and results obtained from the TCL Redheld equation, given 
in Eqs. (lD2|) - flD4j) without the RWA using Eq. (ID5j) and with the RWA using Eq. (1D6[) . 
for three values of the inverse temperature: (a) /3h = 3.0, (b) f3h = 1 . 0 , (c) j3h = 0.5. At 
high temperature, in the QHFP case, the calculations are easier, because there are fewer 
Matsubara frequency terms, while the TCL Redheld calculations are more difficult, because 
more energy eigenstates are needed to account for the high energy excitations. Here we 
included up to 16 states in the TCL case. 


While the QHFP results (red curves) coincide with the exact results (black dots), the 
TCL-Redfield results without the RWA (blue curves) and with the RWA (blue dashed curves) 
are close only near the maximum peak, regardless of temperature. The low-frequency parts 
of the spectra arise from the slow dynamics of the reduced system near the thermal equilib- 
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FIG. 3. Linear response function for the Brownian oscillator system, for three values of the 

system-bath coupling strengths: (a) ( = 0.1, (b) ( = 1.0, (c) ( = 3.0. This function is temperature 
independent in the harmonic case, and we set the inverse temperature to fih = 1. The other 
parameter values are the same as in Fig. [21 The dots represent the analytically calculated exact 
results obtained from Ea. (l50l) . The red, blue, and blue-dashed curves were calculated using the 
real-time QHFP equation, the TCL Redfield equation, and the TCL Redfield equation with the 
RWA, respectively. The intensity of each line is normalized with respect to its maximum peak 
strength. 

rium state, and the discrepancy between the TCL results and exact results arises from the 
equilibration process discussed in Sec. V-A. 


2. Linear response function: Dissipation and non-perturbative effects 

As can be seen from Eq. d^Uji . ipbju;] is temperature independent. Therefore, this function 
is convenient to study the non-perturbative effects of the system-bath coupling, (, and 
non-Markovian effects for slow modulation, controlled by the parameter 7 , apart from the 
temperature effects. In Fig. [3l we compare the linear response functions for the coupling 
strengths (a) ( = 0.1, (b) ( = 1.0, and (c) ( = 3.0 with fixed inverse temperature /3h = 1 
and 7 = 1 . 

While the QHFP results (red curves) coincide with the exact results (black dots), the 
TCL-Redfield results without the RWA (blue curves) and with the RWA (the blue-dashed 
curves) are close only in the weak coupling case considered in Fig. [3](a). For the strong 


23 


















0 1 0 1 0 12 


00 


FIG. 4. The pure non-Markovian effect of Rd)^] investigated in the weak system-bath coupling 
regime. Because the effective coupling strength, ( e ff ~ Cl 2uj o/{l 2 + w o)i depends on 7 , we adjust 
C in each case to keep Ce// equal to its value in the case considered in Fig. [3ja). We chose (a) 
7 = 0.5 and ( = 0.25, (b) 7 = 0.25 and ( = 0.85, and (c) 7 = 0.2 and ( = 1.3 in order to make 
the widths of all the peaks similar. The other parameter values are the same as in the case of Fig. 
[3j The dots, red solid, and blue solid curves are the exact, QHFP and TCL Redfield without the 
RWA results, respectively. 


coupling case considered in Fig. [3] (c), both the QHFP and analytical results exhibit a 
peak near coq = 0.2. This peak arises from the strong coupling between the harmonic 
mode and the low frequency bath mode characterized by 7 2 cu/( 7 2 + cn 2 ) and only appears 
in the simultaneous non-Markovian (7 < cuo) and non-perturbative (£ Wo) case .' 121 The 
existence of this peak, which we call a “non-Markovian bosonic peak, “ is a good indication 
of the applicability of non-perturbative and non-Markovian theories. 

Because the TCL Redfield theory is valid only to second order in the system-bath cou¬ 
pling, the TCL results cannot reproduce this peak. Moreover, the spectrum calculated from 
the TCL Redfield equation without the RWA in the strong coupling case, shown in Fig. 
He), is not positive for cu ~ 5, due to the breakdown of the positivity condition. Despite 
this problem, however, the difference between the TCL results with and without the RWA is 
minor. This is because the spurious behavior caused by the positivity problem is suppressed 
in the non-Markovian treatment of the reduced dynamics, as explained in Appendix B. 
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FIG. 5. The pure non-Markovian effect of investigated in the intermediate system-bath 

coupling regime. We chose (a) 7 = 0.5 and £ = 2.5, (b) 7 = 0.25 and £ = 8.5, and (c) 7 = 0.2 
and £ = 13 in order for the effective coupling strength C, e ff ~ C,l 2(jJ o/{l 2 + w o) t° be the same as 
in the case of Fig. [3Kb) . The dots, red solid, and blue solid curves are the exact, QHFP and TCL 
Redfield without the RWA results, respectively. 

3. Linear response function: Noise correlation and non-Markovian effects 

We next discuss the non-Markovian effects in the Brownian oscillator system. It should 
be noted that when the inverse noise correlation time, 7 , is decreased, the effective coupling 
strength becomes stronger, even if we fix £, because the bath can interact with the system 
multiple times when the correlation time is long. In order to study the pure non-Markovian 
effects, here we employ an effective coupling strength £ e ff ~ <5r 2 (u;o) = £ 7 2 cuo/( 7 2 + cUq)^ 
and fix it while varying 7 . 

In Fig. [0 we plot in the weak coupling regime corresponding to Fig. [3Ka) . Here¬ 

after, we do not consider the TCL Redfield equation with the RWA, because the difference 
between the TCL results with and without the RWA is minor. While all of the peak profiles 
are similar if we fix £ e //, the peak position shifts slightly in the high-frequency direction, 
because a change of 7 results in a change of 5h2 2 (a;o). As the exact results and the QHFP 
results in Fig. [4] indicate, there is no clear indication of non-Markovian dynamics in this 
weak coupling regime, once we have normalized the effective coupling strength. 

While the TCL Redfield results are close to the exact results in the fast modulation (weak 
non-Markovian) case depicted in Fig. [4](a), they differ significantly in the slow modulation 
(strong non-Markovian) case considered in Fig. [4[c). This is because the perturbative 


25 


















3 



<»! 


( 2 ) 

FIG. 6. The second-order response function R^j, r [uj 1 ,^ 2 ] of the Brownian oscillator system cor¬ 
responding to the intermediate coupling case considered in Fig. 0 (b). The results here were 
obtained from (a) the analytical expression Ea. (|54p . (b) the QHFP approach, and (c) the TCL 
Redfield approach without the RWA. The intensity of each graph is normalized with respect to the 
maximum peak intensity. 


description of the TCL Redfield equation breaks down as a result of the fact that multiple 
system-bath interactions arise due to the slow modulation, even in the weak coupling case. 
Thus the TCL-Redfield result without the RWA becomes negative for oj > 4. 


In Fig. 0 we plot j n the intermediate coupling regime corresponding to Fig. [3](b). 

It is seen that while the QHFP results always coincide with the exact results, the discrepancy 
between the TCL Redfield and exact results is large in the slow modulation (strong non- 
Markovian) case, due to the non-perturbative nature of the interactions. Specifically, the 
lack of a non-Markovian bosonic peak becomes apparent even at this intermediate coupling 
strength if the modulation is slow. Moreover, the TCL result without the RWA becomes 
negative in the region uj > 2.2. Because the non-Markovian effects in dynamics make the 
non-perturbative nature of the interaction conspicuous in the case of slower modulation, 
due to the existence of the multiple system-bath interactions for slow modulation, the TCL 
Redfield equation does not have the capability of treating pure non-Markovian effects even 
it reproduces the high-frequency part reasonably well. 
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C. Nonlinear response function: Dynamical system-bath coherence 

As explained in Appendix C, the system-bath interaction induces static effects arising 
in imaginary time and dynamic effects arising in real time and complex time. While the 
static effects can be obtained from the thermal equilibrium distribution, as illustrated in Sec. 
V-A, we have to study the nonlinear response function to elucidate the dynamic effects. It 
should be noted that, in addition to their inability to treat strongly non-Markovian dynam¬ 
ics, the conventional reduced equation of motion approaches involving the TCL Redfield 
equation have a severe limitation in studying systems subject to time-dependent external 
forces because their description of the damping kernels is based on energy eigenstates.— The 
capability of an approach to treat external forces can also be examined by calculating nonlin¬ 
ear response functions, because nonlinear measurements can capture the effects of multiple 
interactions through time-dependent external forces. Here, we calculate the second-order 
nonlinear response function of the position given by 



(53) 


This is an observable in two-dimensional THz-Raman spectroscopy system. 1221123 Note 
that, because of the Gaussian integral involved in the expectation value ((•••) = tr{- ■ 


■ exp(— /3H to t)}), the contribution from the lowest-order response, ([[g(fi + t 2 ), <?(U)], <?])/h 2 , 
vanishes.- 1 ^ In the harmonic case, there is also a contribution from R^ RT (ti,t 2 ) = 


— ([[g(ti + t 2 ) ,g 2 (ti)] ,q])/k 2 , which corresponds to an observable in 2D THz-Raman-THz 
spectroscopy system. We find that to explore the system-bath coherence, Eq. (l53l) is suitable, 
as we show below. This response function in the harmonic Brownian case can be calculated 
analytically as^ 


-^ ttr (^ i , ^ 2 ) — C(t 2 )C(ti + t 2 ), 


(54) 


where C(t ) is obtained from the Fourier transform of Eq. d49lh To apply the Liouvillc operator 
formalism, we rewrite Eq. (j53l) as 



(55) 


Using the above expression, we calculated Rx RR (ti,t 2 ) for various values of t\ and t 2 by 
extending the method employing Eqs. (l5Tf and (l52]b &^£ The response functions evaluated 
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from Eqs. (1541) and fl55]) are then Fourier transformed to obtain two-dimensional spectra, 
Rttr [kb j ^2] • 

In Fig. El we plot 2D spectra in the frequency domain obtained from (a) the analyti¬ 
cally exact approach, (b) the QHFP approach, and (c) the TCL Redfield without the RWA 
approach under the same physical conditions as in Fig. [5] (b). We find that while the an¬ 
alytically exact and QHFP results exhibit peaks at (a^i,a^ 2 ) = (0,1) and (tni, u; 2 ) = (1,1), 
the TCL approach cannot reproduce them. As shown in a study of multi-dimensional 
spectroscopy, in order to have these peaks, the dynamical system-bath coherence sub¬ 
ject to the second interaction at time t\ must be maintained throughout the time evolu¬ 
tion described by Eq. (|55]h — I 11 the TCL case, however, the time evolution is described in 
terms of the reduced operator tr B {G(t)}, derived from the factorization assumption with 
tr ^q 2 tr B {G(t2)}q x tr B {G(ti)}q x tr B {p^ t }^. While the exact dynamics maintain the coher¬ 
ence during the period of length t\ + £2 expressed by C{t\ + £ 2 ), the TCL approach cannot 
maintain this coherence. In contrast to the Redfield approach, because the HEOM approach 
can store this coherence in the hierarchal members, it is capable of treating a nonlinear re¬ 
sponse function. 

Because many modern experiments utilize the nonlinear response of a system, which is 
measured by applying a variety of time-dependent external forces, the capability to cal¬ 
culate the nonlinear response function is important. The validity of the HEOM approach 
has been demonstrated for systems subject to time-dependent external forces — ~ 97 4 05 In 
addition to the HEOM approach, the path integral approach has also been shown to have 
this capability.— 


D. Thermal equilibrium state and thermodynamic quantities 

We finally examine the imaginary QHFP equation by considering our results obtained 
through numerical integration of Eq. fj36j) from r = 0 to / 3h using the harmonic potential 
to compare W^ g Q (p,q) presented in V-A and the partition function Z, 4 . The number of 
Matsubara frequencies used in the imaginary-time QHFP is K = 4. The mesh size was 
optimized for the Euclidean Liouvillian, and we chose n p = 60-120 and n q = 120-240. 
Because the distribution is spread relatively widely in the higher temperature case, we 
employed a coarser mesh in that case. 



(a) z —0.0 


(b) t —0.1 



FIG. 7. Solution of the imaginary HEOM at four values of the imaginary time, r. Here, we plot the 
zeroth member, q ; r), only. The initial state is presented in (a) r = 0, while the final state 

is presented in (d). We confirmed that the normalized distribution of the state in (d) is identical 
to the distribution given by Eci. (|42|) . within numerical error. 



P 


FIG. 8. The partition function, Za , entropy, Sa , internal energy, Ua, and heat capacity, Ca , of a 
Brownian oscillator system calculated using the imaginary-time QHFP as functions of the inverse 
temperature, /3h. The dotted curve represent the partition function obtained from the analytical 
expression, Ea. (f56l) . Because analytically calculated exact results and the HEOM results, Za, are 
nearly identical, here we plot only Sa, Ua, and Ca for the HEOM case. 
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In Figd we display solution of the imaginary-time QHFP, Eq. (j36j) with (3h = 1 for several 
values of r. Because the damping kernels in the imaginary-time QHFP are defined by the 
Matsubara frequency at f3h , the solutions r < /3h do not correspond to the equilibrium dis¬ 
tribution at temperature r. While the initial distribution is flat, the distribution approaches 
a Gaussian form due to the Euclidean and the damping operators. At t = /3h, the solution 
coincides with the analytical solution given in Eq. (l42]l with Eqs. (l45j) and (l46]l . 

While the equilibrium distribution can also be obtained from the real-time QHFP, 
as shown in Sec. V-A, the thermodynamic quantities can only be calculated from the 
imaginary-time QHFP. We next demonstrate this point. In the BO case, the partition 
function can also be evaluated analytically in terms of the Matsubara frequencies as^- 

-i 00 2 

= n 2 1 2 U 1^2( Y ( 56 ) 

f3ftw 0 ^ + + 8T 2 (u k ) 

We should note that, the normalization constant of the real-time QHFP is N = 2n^/ ( p 2 )(q 2 ), 
whereas that of the imaginary-time QHFP is Za obtained from Eq. (l40l) . Because Zbo 
involves a temperature dependent factor other than N, we cannot calculate the partition 
function using the real-time HEOM approach. 

To obtain thermodynamic quantities, we first repeated the integration of the imaginary- 
time QHFP from j3h = 0.025 to 3.05 with step size A/3h = 0.025 to derive Za- Then, 
we calculated thermodynamic quantities through Za- In Fig. [8j we compare the partition 
function given by Eq. (I40l) (brown curve) and that obtained from Eq. (1561) (dotted curve). 
As expected, the imaginary-time QHFP results coincide with the exact results. For the 
purpose of demonstration, we also plot the entropy, Sa, the internal energy, Ua, and the 
heat capacity, Ca, calculated with the imaginary-time QHFP. The behavior in the high 
temperature regime is very different from that in the spin-Boson case,— because the BO 
model has an infinite number of excited states. 

VI. CONCLUDING REMARKS 

In this paper, we presented real-time and imaginary-time QHFP equations derived using 
the influence functional formalism with correlated initial conditions. While we found that 
the QHFP equations in real time possess the same form as those obtained from a factorized 
initial state, we introduced a modified terminator in order to facilitate the more efficient 
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calculations of non-Markovian dynamics. 

The capability of the real-time QHFP was verified through non-perturbative and non- 
Markovian tests based on (i) the steady-state distribution, (ii) the symmetric auto-correlation 
function, (iii) the linear response function, and (iv) the nonlinear response function. This 
was done to test the capability of the real-time QHFP to properly model the effects of (i) 
static system-bath coherence, (ii) fluctuation, (iii) dissipation and non-Markovian effects, 
and (iv) dynamical system-bath coherence, respectively. The ability of the model to ac¬ 
count for the dynamical system-bath coherence is particularly important if we wish to study 
dynamics under time dependent external forces. While many of the methodologies devel¬ 
oped for reduced quantum dynamics have been tested only with regards to the relaxation 
dynamics of the population state over short periods of time, the long-time behavior of the 
dynamics, represented by the low frequency parts of the correlation functions, is essential 
to test the capability of this approach for non-Markovian dynamics. Because the bath can 
interact with the system many times in the case of slow modulation, the dynamics of the 
reduced system can only be described with a non-perturbative treatment when the system is 
strongly in non-Markovian. For this reason, the non-perturbative treatment and the mixed 
state (or unfactorized) treatment of the system-bath interactions are both important. 

In this paper, we considered only the harmonic case, the HEOM approach can be used 
to treat potentials of any form with time-dependent external forces. Although it had not 
been shown until the present paper that the QHFP equations derived from correlated ini¬ 
tial conditions have the same form as those obtained from factorized initial conditions, the 
usefulness of the real-time QHFP approach has been demonstrated for various problems in¬ 
volving chemical reactions,^’ 95 photo-dissociation, 96,9 ' nonlinear optical response ,— 8 ~— 02 res¬ 
onant tunneling , 103,104 quantum ratchets,— and tightly bound electron-phonon system.— 
However, with the modified terminator introduced in this paper, the same calculations can 
be carried out more efficiently. 

A confined potential system involving a Brownian oscillator system can also be treated 
using the HEOM approach in the energy eigenstate represent at ionA2£ in the same manner 
as in the present study of the TCL Redfield equation, but quantum transport problems 
characterized by open or periodic boundary conditions can be studied only with the QHFP 
approach,—because we cannot introduce the energy eigenstates for this kind of 
problem. Nonlinear system-bath coupling, which plays an important role in vibrational 
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spectroscopy, can also be taken into account in the QHFP formalism.—— - -^ Extension to 
multi-potential surfaces is also possible.—' 9 ' Because this formalism treats the quantum 
and classical systems with any form of potential from the same point of view, it allows 
identification of purely quantum mechanical effects through comparison of classical and 
quantum results in the Wigner distribution.—' 102 ! 105 

We showed that the thermal equilibrium state obtained from the imaginary-time QHFP 
is equivalent to the steady state solution of the real-time QHFP. Because the imaginary-time 
QHFP is defined in terms of integrals carried out over the definite time interval, we were able 
to calculate the equilibrium state more easily in this case than in the case of the real-time 
QHFP. Moreover, using the imaginary-time QHFP, we were able to calculate the partition 
function, and from this, we could directly obtain several thermodynamic quantities, namely, 
the free energy, entropy, internal energy, and heat capacity of the system in the dissipative 
environment. Numerical integration of the real-time and imaginary-time QHFP equations 
is computationally intensive. Nevertheless, we were able to study the dynamics of one¬ 
dimensional potential systems using personal computers.— - — Great effort has been made to 
reduce the computational intensiveness of algorithms used to implement the real-time HEOM 
approach. For example, the hierarchy has been optimized for numerical calculations,— - — 
and a graphic processing unit (GPU) 132 and parallel computers^ have been utilized in 
order to facilitate the treatment of larger systems and to treat non-Drude type spectral 
distribution functions.— - — The same techniques can be applied to the case of real-time and 
imaginary-time QHFP equations. 

As supplementary materials, we supply the FORTRAN codes for the real-time and 
imaginary-time QHFP, entitled TanimuranFP15 and ImTanimuranFP15, to help further 
development in this field.- 131 
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Appendix A: Influence functional with correlated initial conditions 


The reduced density matrix elements of the system are obtained in path integral form as 




J B 


rq=q(t) 


J tot J q Q =q(0) 


D[q(t )] 


rq'=q'{t) 


lq' 0 =q'{0) 


D[q'(f)} 


i-q^qiPH) 


lq o =q(0) 


D[q(T )] 


X e ji s A c tA e ®[<h<i', c i; t,ph] e -^S A [q;T] e -^S A [q',t] 


(Al) 


where Zb is the partition function of the bath and the Euclidean action is given by 

Sa[<!;t\ = J dr !^mq{r) 2 + U(q(r))\. (A2) 

The influence functional for the correlated initial state expressed in terms of the influence 
phase is given by— 


$[?,?', 9 ; t,pn] = 


; \ 2 ,t 


,ih 


dt"^BpS)q x {t")q 0 {t"^ 


h 

• \ 2 rt 


0 


- (~j) J Q dt "J o dt , q x (t , ')[-tL 1 (t ,, -t , )q%t , ) + L 2 (t ,, -t , )q x (t'^ 


+ 


k 2 , 
1 


ft rph 

/ dt" dr'q x (t")q{r')L(t" + it') 

o Jo 

r/3fr i /"PS 


dr" / dr'q{r")q(r')L (r" - r'), 


(A3) 


where q x (t) = q{t) — q'(t ) and q°(t ) = q(t ) + q'{t). Using the spectral density, we 

rewrite these functions for 0 < r < /3h as 


2 

L(t + fr) = — / dujJ(uj ) 

Ph Jo 


1 f 2w 

- + -.2 , ..9 cos(z/ fc r) 


CU ' l>l + U 2 
k =1 


cos (cut) 


(A4) 


, . 2 f°° , „ 2 

+*— / duJ(uj) > „ , „ 

PKJ o + 


In the case r = 0, we have L(t) = iL\{t) + L 2 (t) with 

poo 

Li(t) — / dujJ{u) sin(a;f), 


sin(z/ fc r) sin(u;f). 


(A5) 


L 2 (t ) 


dujJ(oj) cotli 



cos(u;f), 


(A6) 
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and in the case t — 0, we have L{r) = L{ir) with 

2 f°° 

L(T) = Jhl 

where the quantities V}. = 27T k//3h are the Matsubara frequencies. For later convenience, we 
also introduce the canonical correlation 


1 

id 


E 

k =1 


2u 


vt 


■ 


cos (u k r) 


(AT) 


, . 2 f°° , J(w) . . 

B(t) = — / du --cos(a jt), 

ri Jo id 


(A8) 


and express the counter-term of the potential using 5(0). 

The function L 2 (t) is related to B(t) through the quantum version of the fluctuation- 
dissipation theorem, L 2 [uj\ = Hu coth(/3fruj/2)B[uj]/2, which insures that the system evolves 
toward the thermal equilibrium state, trB{exp[— {3H to t}}, for finite temperatures in the case 
that there is no driving force 
Using the relations 


- / dt"B(0)q x (t")q°(t") + / dt 


, dB(t" — t') 


i dt r - V (t")g°(0 


dt' 


rt" 


= - / dt"B(t")q x (t")q°( 0)- / dt" / dt' B(t" -t')q x {t") 


dq°(f) 

dt' 


(A9) 


and 


1 f 0h 

»L dT " 


1 f 0h 

dr'q(r")q(r')L (r" - T ') -—J dr"B(0)f(r") 


1 

2h? 


r/3h 


'•ph 


dr" / dr'q(r")q(T')L (r" - r '), 


(A10) 


the influence functional can be rewritten as 


F C i[q,q',q-, t,ph] 



x e( 
x e( _ 
x e* 1 


ji) 2 fo dt"^B(t")q x (t")q°(0) 
if f* dt"qX (t") if dt’ f B(t"-t’) a -^gl 
l) 2 /o (*") if dt’L 2 (t"-t') q x (f) 

tidt"CdT’qX(t")q(T')L(t"+iT') 


(AH) 


where we have included the bath part in the initial thermal state of the system p e ^[q] (3h] as 
p eq [q; ph] = Z B e-* SA[q;T]+ ^ J ‘ dT ''C dT '^ T ")^W-r')_ ( Al2 ) 
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FIG. 9. The roles of the system-bath interactions illustrated schematically. The solid black lines 
repr esent the wave function of the system A along the complex counter-path (see Fig. 1 in Ref. 


108l ), and the arcs and blue lines represent the system-bath interactions. Because the bath degrees 


of freedom have been reduced, the bath interactions connect the wave function of the system A 
at multiple complex times. The blue arcs and lines correspond to the fluctuation and dissipation 
processes described by the terms containing B(t" — t') and L 2 {t" — t') in Ea. dAllD . while the red 
arcs represent the static thermal system-bath correlation described by L{t" — t') in Ea. (lA12D . The 
green arcs represent the correlation in complex time described by L(t' + t'), which leads to the 
correlated initial conditions. 


The contributions arising from the factorized initial conditions or correlated initial conditions 
consist of two parts. One is a static contribution represented by the term containing the 
imaginary-time integrals of L (t" — t') in Eq. (jAl2jl . Because of this term, the thermal 
equilibrium state of the system is not the equilibrium state of the system alone (pure state), 
but that of the combination of the system and bath (mixed state). The other is the correlated 
state contribution, represented by the term containing the complex time integrals of L(t' + 
it') in Eq. flAlljh The second contribution involves the effects of the dynamical correlation 
and is negligible when the Markovian assumption is applied, while the first contribution 
always plays a significant role. It is important to note that, in addition to the fluctuation 
and dissipation denoted by L 2 (t ) and Bit), respectively, there is a dynamical contribution 
from the correlated initial conditions. The role of the system-bath interaction is illustrated 
in terms of each contribution in Fig. [9j 
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FIG. 10. The noise correlation function, L 2 (t), depicted as a function of the dimensionless time 
t for several values of the inverse noise correlation time: (a) 7 = 0.5, (b) 7 = l,(c) 7 = 5 . Note 
that 7 —> 00 corresponds to the Markovian (Ohmic) limit, as can be seen from Eq (|5j). The inverse 
temperatures are, from top to bottom, f3h = 0.5, 1.0, 3.0, and 5. The noise correlation becomes 
negative in (b) and (c) at low temperature (large f5h ) due to the contribution of the Matsubara 
frequency terms. 

Appendix B: Drude spectral distribution and the violation of the positivity 
condition 

With Eq.([5]), for 0 < r < /3h, we obtain^ 

{ OO 

c" + Wk cos (u k r) + ic' k sin(z/ fc T)] 

k= 1 

OO 

+ ^2 c ' k [ cos ( z/ fc r ) _ i s^M] e ~ Ukt i 

k =1 

where c'q = mCl/P, c' k = — 2m(^ 2 u k /— v k ), and c k = 2m('y 3 //3(7 2 — v k ) for k < 0. At 
t — 0, the above equation reduces to 

OO 

L( t ) = y^c fc cos(t/ fc r), 

k =0 

where u 0 = 0, c 0 = Cq, and c k = c' k + c" for 1 < k, while at r = 0, we have 

B(t) = m( r ye~" /t (B3) 


(B2) 
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and 


OO 


MS) = + £V t e--*‘. (B4) 

k =1 

As shown in Fig. [TUI the noise correlation, L 2 (t), becomes negative at low temperature. 
This results from the contribution of the terms with u k = 2nk//3h in the region of small t. 
This behavior is characteristic of quantum noise. 8 We note that the characteristic time scale 
over which we have L 2 (t) < 0 is determined by the temperature and is not influenced by the 
spectral distribution J(oj). Thus, the validity of the Markovian (or <5(t)-correlated) noise 
assumption is limited in the quantum case to the high temperature regime. Approaches 
employing the Markovian master equation and the Redfield equation, which are usually ap¬ 
plied to systems possessing discretized energy states, ignore or simplify such non-Markovian 
contributions of the fluctuation, and this is the reason that the positivity condition of the 
population states is broken.— — 

As a method to resolve this problem, the rotating wave approximation (RWA) (also 
known as the ’’secular approximation”) is often employed, but a system treated under this 
approximation will not satisfy the fluctuation-dissipation theorem, and thus the use of such 
an approximation may introduce significant error in the thermal equilibrium state and in 
the time evolution of the system toward equilibrium. Because the origin of the positivity 
problem lies in the unphysical Markovian assumption for the fluctuation term, the situation 
is better in the non-Markovian case, even within the framework of the Redfield equation 
without the RWA, as discussed in Sec. V. In the classical limit, with h tending to zero, L 2 {t ) 
is always positive. 

While conventional approaches employing reduced equations of motion eliminate the 
bath degrees of freedom completely, the HEOM approach retains information with regard 
to the system-bath coherence in the hierarchy elements. Because of this feature, the HEOM 
approach can treat the reduced dynamics in a non-perturbative, non-Markovian manner. 
To obtain a more compact form for the HEOM, we use the following approximate form 
for L 2 (t), given in Eq. (JB4J) : L 2 (t) ~ + Ylk=i cf k e ~ Vht + E^a'+i with 

Cq = hm^ 2 cot(/3^7/2)/2. Here, we choose K so as to satisfy v k = 2nK/(/3h) 3> cu c , where 
ui c represents the characteristic frequency of the system. Under this condition, we can apply 
the approximation v k eT Vk ' t \ — C k 5(t) (for k > K + 1) with negligible error at the desired 
temperature, 1//3, where C k = ^k/( u k + w c) the correction factor that compensates for the 
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overestimation of L 2 (t) in the approximation at very low temperature for small cut-off K. 
The accuracy of this approximation is verified on basis of the asymptotic behavior of L 2 (t) 
as a function of K. Then, the HEOM can be obtained by considering the time derivative 
of Eq. (J5]).~ When the temperature becomes high (i.e. for j3ftr/ <C 1), the noise correlation 
function reduces to L 2 (t ) ~ and hence the noise modulates the system as a 

Gaussian-Markovian stochastic process.— 1 ^ 


Appendix C: Derivation of the HEOM in configuration space 

In the present appendix, we construct the equation of motion for Pjl] 2 ...j K (q, ( l'\ t). In order 
to obtain differential equations in time, we consider the reduced density matrix elements at 
t + St, 


. . i r r rq(t)=q-y rq' (*)=?' -y' 

phh-jK^^'’ t + 6t ) = 42 / d v / d v‘ £%(*)] / D Wif)} 

J J J q(0)=qo Jq'(0)=qo' 

rq' 0 =q(PH) 

x/ D[q(T)]p*[i,ffi 

Jqo=q(o) 

t-\- 5 t -i 

dt! e* 7© 0 (t') + G 0 (0) - -S(l3h) 

L</0 h 

t+St i 1 7 3k 


x |e 

K 

x n <' 

k =1 


-7 (t+5t) 

-v k (t+8t) 


L</0 


n 


X e iS[q,t+5t] Fci ^ ^ y. t + 5tj p h } e -iS[ g ',t+6t]' 


(Cl) 


where A is the normalization constant for the integrals over y and y', and we set q = q(t ) + y 
and q' = q'(t) + y' with q(t + St) = q and q'(t + St) = q'. We then expand Pj^j 2 ...j K (Q, Q r ] t + St) 
in terms of 5t up to first order. Because y and y' also depend on St, we have to expand the 
above equations in terms of y and y'. In the following, we expand the components separately. 

The action part can be expressed as 

e iS{q,t+8t] = A[f(ft) 2 - U (l-y')] St e iS[q-y,t] 


= e* ( 1 - 'A-U (,) ) 


(C2) 
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The influence functional is evaluated as 


F C i[q,q',q-, t + 6t,pK} = 


1 - St $(t)\ e“ 7( * +ft) 


pt+8t 


Uo 


dt' e 1 *' ^/@ 0 (t') + G 0 (0) — —Q(/3h) 

n 


K 


~ St $(t) l J^e 


-u k (t+St) 


k =1 


L/0 


t+8t -j 

dt’e^' v k Q k {t')--^ k {ph) 


x F C i[q -y,q' - y',q; t,/3h]. 

In the following, we apply the Gaussian integrals 

1 / imy ^ 

- / dyye^nsr = o 


and 


1 f o im 2 Xfi 

- / dyy eawt 1 ' = — St, 
A m 


- StZ'(t ) 
(C3) 

(C4) 

(C5) 


where the normalization constant is chosen to be A = f dy exp(imy 2 /2h5t). Gaussian 
integrals higher than fourth order can be ignored, because they produce contributions smaller 
than o(5t). 

With Eqs. (1C2|) and fjC3|l . the expansion of the last term in Eq. (1CID is completed by the 
following: 

= (l - 4 - 4 + + V Tw) 

x ei s ^F CI [q, q', q\ t, /3/f]e4 5 (<?'>*). (C6) 

Then, collecting the pieces from Eqs. (1C2|) . ()C3[) and (1C6[) . and keeping terms up to o(St), 
we have the following for the kinetic term of the Hamiltonian: 


dy imy 2 f d y 2 d 2 
e 2Mt ( 1 — y ——|- 


A 


= 1 - St-. 


d_ 2 _ 

h V 2 m dq 2 


(C7) 


dq 2 dq 2 / 

We next consider the expansion of the factor {• • • } n in Eq. (1C 1 D . first in terms of y and 
y', and then in terms of St. For the expansion in y and y' up to second-order, we have 


(y + y') { f d^'7e- 7(^ - ^,) 0o(^ , ) + Go(0) - l&m 


h 


n—1 


(C8) 


This term reduces to (q, q';t), and therefore the contribution of the above to the 

relevant order in St can be expressed as 


nmC'y f dy f dy' im y 2 imy' 2 ( d . d y 2 d 2 y' 2 d 2 

I / -x e2Mte " 2Mt \ 1 ~ y ^- y ^ + 


dq dq' 2 dq 2 2 dq' 2 

x(y + y') (<?, q\ t) 


(C9) 
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Then, integrating over y and y', we have 


-St 


inh( , y 


d d 


PhJL 


2 [ dq dq 1 

For the expansion of {•••}"" in terms of St, we have 

n0 o (f){ fc^e-^-^eotO + G 0 (0) - 


(CIO) 


71—1 


St 


-ri7i / dt'-fe-^-^Qoit') + G 0 (0) - -O(Ph) \ St. 


(Cll) 


We can expand the factors {• • • p fc in Eq. (lCl|) similarly to the {• • • } n factor. We obtain 

1 


/ t 1 'j 3k 

dfVfce-^O-CQ^') _ _«|r k {j3h)\ St 

-J^ k e k (t)^-fyt'u k e-^-ne k ^ _ | J St 

Using the definition of the hierarchy elements Eqs. (fl4l) and we obtain 

{-^7(<?, q'] t ) + nQ 0 (t)pf~ l ] K (q, </; f)} St 

and 


(C12) 


(C13) 


-jk^kP^l.. dK (q, gt) - ]kip^hi>)i>y'L. j , i . jK (q, q P } St, (C14) 


From Eqs. (lCll|) and (IC12|) . respectively. Finally, substituting the results from each of the 
above expansions, contained in Eqs. (lC7D . (1C13D and (jC14f) . into Eq. (lCl|) . we construct the 
complete form for this expression to o(St), and from this, we obtain Eq. fflGl) . 


Appendix D: Time-Convolutionless (TCL) Redfield Equation 

The TCL Redfield equation is the reduced equation of motion in the case of non- 
Markovian noise whose damping kernels are expressed in a time-convolutionless form.— 
The TCL Redfield equation is exact if the system Hamiltonian, Ha, is time-independent 
and if Ha commutes with the bath interaction. However, for the BO model considered in 
this paper, defined by Eq. (JTJ) , the system Hamiltonian does not commute with the bath 
interaction. 

In order to apply the Redfield theory, we need to use the eigenstate representation of 
the system. For this reason, we include a counter-term in the system Hamiltonian, and 
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consider the modified Hamiltonian, H' A = {H A + m( , yq 2 /2). For the jth eigenenergy, E'- = 
h(j + 1/2 )ojq, where u' 0 = \/oJq + CTj the eigenfunction for H A is expressed in terms of 
Hcrnrite polynomials, Hj(-), as 


^M) 






(Dl) 


where a = sJrruJjh. We denote the ket vector for ipj(q) by | j). The TCL Redfield equation 
for the reduced density matrix elements, Pjk(t) = (J\pA(t)\k), is then given by 

d 

~Q^Pjk(t) ^^jkPjk{t) ~\~ ^ ' Rjk,lm{t) Plmjt), (D2) 

l,m 


where E jk = (j — k)u>' 0 and Rjk,im{t ) is the Redfield tensor defined by 


Rjk,lm(t ) — r mk,jl(t) + rj j Sjm b j n ,nl(t) $jl ^\n,nrrSf)i 


(D3) 


with 


1 - e -(7+Mm) t 


Tjfc,Zm(^) Tj'fcZm, 


— y 

Rh ^ 


C7 2 rv 1 e -K'+^Lh' 


2 sin(®^) 1 + Hm Vk'+Mm 


• (D4) 


The interaction tensor is defined by 


Tjkim = (j\q\k)(l\q\m). (D5) 

The rotating wave approximation (RWA) is expressed as q£j = \/2h/muj' 0 (a + + a~)(bj + 
bj~) ~ d + bj "b where a ± and b± are the creation and annihilation operators of the BO 

oscillator and the jth bath oscillator, respectively. For (j\a + \k) = 0 and (k\a~\j) = 0, with 
j ^ k — j + 1, the interaction tensor in RWA form is given by 

9b 

r furn = ~~F ({j\a + \k){l\ar\m) + ( j\a~\k){l\a + \m )) . (D6) 

TTIUJq 

In the Wigner representation, the eigenstate elements of the density matrix are expressed 
as 

W jk {p, q) = ^ J dx cos (y) ipj (q - Va (? + |) • (b>7) 

The total distribution is then given by 

M 

W(p,q\t) = P jk {t)W jk (p,q), (D8) 

j,k =1 
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where M is the number of energy eigenstates employed to solve the TCL Redheld equation. 
The factorized initial state is expressed as 

M 

W(p, q; 0) = V — exp(-f3E / j )W n (p, q), (D9) 

j=i 

where Z' A = . exp (—j3Ej). By comparing the steady state solution of the TCL Redheld 
equation in the Wigner representation with the analytical solution of the BO model given by 
W^q (p, q) with Eqs. (145|) and (l46]h we can check the accuracy of the steady-state distribution 
in the TCL formalism. 
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